% Table for the missing city
% 2017-2-27
% Tentative plan: We provide figures for houston
clc; clear all; close all;
workpath00 = pwd;

% predX_M1: covariates at the tract centroid is also estimated via
% hierachical model 
% M1: includes intercept and log(1+dist), coefficients are varying with Area
% M2: includes intercept and log(1+dist) and log(dist 2 coast), coefficients are varying with Area
estXmode = 3; %1=M1, 2=M2, 3=M3


%% Figure option
cityid = 520; %520 for Atlanta, 1600 for Chicago, 4480 LA, 5600 NY, Houston 3360
nyear = 6;
modelind = 6; %2: linear, 6: quartic, size-cubic
islog  = 2; %2 for log(1+x)

%% Load Data and preliminary stuff
load('data_all');
% load('results_full_logp1.mat');
% load('results_full_logp1_A.mat');
% load('results_full_logp1_A_s3_REML.mat');

% % other info
% ncity = size(ZZ2,1);
%
% % Spatial function
% ngrid = 100;
% xgrid = linspace(0,60,ngrid)';
% logxgrid = log(1+xgrid);
%
% logxgrid1 = log(1+xgrid);
% logxgrid4 = [logxgrid1, logxgrid1.^2, logxgrid1.^3, logxgrid1.^4];

%% Missing cities
% [120; 1350; 1480; 1580; 1620; 3400; 3850];
pmsaname_ms = {'Albany, GA', 'Casper, WY', 'Charleston, WV', 'Cheyenne, WY', ...
'Chico-Paradise, CA', 'Huntington-Ashland, WV-KY-OH', 'Kokomo, IN'};

% 120 Albany, GA, cityhall: 31.577734, -84.152437
% 1350 Casper, WY cityhall: 42.851573, -106.327282
% 1480 Charleston, WV cityhall: 38.355842, -81.639657
% 1580 Cheyenne, WY city hall: 41.135669, -104.822013
% 1620 Chico-Paradise, CA, city hall: 39.728895, -121.838067
% city hall 1: 39.729147, -121.837649
% city hall 2: 39.749591, -121.634307
% 3400 Huntington-Ashland, WV-KY-OH
% city-hall 1: 38.419628, -82.444956
% city-hall 2: 38.478755, -82.637708
% 3850 Kokomo, IN: city hall 40.486118, -86.129631
% load updated data
% save('censustracts4_ch.mat');

% load('censustracts4_ch.mat');
% tract_ms.combofips = combofips;
% tract_ms.ch_lat = ch_lat;
% tract_ms.ch_lon = ch_lon;
% 
% tract_ms.lati = lati;
% tract_ms.loni = loni;
% 
% tract_ms.min_dist = min_dist;
% tract_ms.landsqmi = landsqmi;
% tract_ms.landsqmi_ua = landsqmi_ua;
% tract_ms.d2c = dist_tract_coast;
% save('tract_ms.mat', 'tract_ms');


load('tract_ms.mat');

%% Load esitmation results
load('results_full_logp1_A_s3_REML.mat');
est_m = est_s6; %final model
temp_est = est_m;

load('results_estX_REML.mat'); %model for covariates
if estXmode == 1
    predX = estX_M1.predX;
    predM = estX_M1;
elseif estXmode == 2
    predX = estX_M2.predX;
    predM = estX_M2;
elseif estXmode == 3
    predX = estX_M3.predX;
    predM = estX_M3;
end



%% Compute land values for missing cities

xtflag_ms = tract_ms.combofips;
xtval_ms = [120; 1350; 1480; 1580; 1620; 3400; 3850];
ncity_ms = length(xtval_ms);


% aggregate info
temp_x = XX;
temp_y = YY;


m_XX = mean(XX);
m_XX2 = mean([XX(:,14).^2, XX(:,14).^3]); %nonlinear term in size
npol = size(temp_est.post.m,1)-6; %number of polynomials


temp_min_dist= log(1+ tract_ms.min_dist);
grid = [temp_min_dist, temp_min_dist.^2, temp_min_dist.^3, temp_min_dist.^4];
grid = grid(:,1:npol);

% loop start here
% -------------------------------------

%
tab_area = [];
tab_ua   = [];
tab_land = [];

tab_lv = [];
tab_name = [];

for j = 1:1:ncity_ms;
    temp_combofips = xtval_ms(j);
    
    
    % relevant information
    temp_sel = (tract_ms.combofips == temp_combofips);
    temp_n = sum(temp_sel);
    temp_d2c = tract_ms.d2c(temp_sel);
    temp_ua = tract_ms.landsqmi_ua(temp_sel);
    temp_ZZ2 = [1, log(sum(temp_ua))];
    
    % construct prior for missing city
    b0 = [...
        est_m.a0, est_m.a1; ...
        est_m.td(1), 0;
        est_m.td(2), 0;
        est_m.td(3), 0;
        est_m.td(4), 0;
        est_m.td(5), 0;
        est_m.d10, est_m.d11; ...
        est_m.d20, est_m.d21; ...
        est_m.d30, est_m.d31; ...
        est_m.d40, est_m.d41; ...
        ];
    
    % prior mean and covariance
    m0 = b0*temp_ZZ2';
    V0 = est_s6.post.Sig2;
    
    
    % temporary variables
    cityid = xtval_ms(j);
    
    temp_sel  = tract_ms.combofips==cityid;
    temp_grid = grid(temp_sel,:);
    temp_lon = tract_ms.loni(temp_sel);
    temp_lat = tract_ms.lati(temp_sel);
    temp_d2c = tract_ms.d2c(temp_sel);
    
    temp_area = tract_ms.landsqmi_ua(temp_sel); %urban area
    temp_area_n = temp_area/sum(temp_area); % normalized area
    
    temp_area2 = tract_ms.landsqmi(temp_sel); %urban area
    temp_area2_n = temp_area2/sum(temp_area2); % normalized area
    
    temp_n = size(temp_grid,1);
    
    % constructing predicted covariates
    temp_pred_ZZ = repmat(temp_ZZ2(1,2), temp_n, 1);
    if estXmode == 1
        new_x = [ones(temp_n,1), temp_pred_ZZ, ...
            temp_grid(:,1), temp_pred_ZZ.*temp_grid(:,1)];
        new_z = {[ones(temp_n,1), temp_grid(:,1)]};
    elseif estXmode ==2
        new_x = [ones(temp_n,1), temp_pred_ZZ, ...
            temp_grid(:,1), temp_pred_ZZ.*temp_grid(:,1), ...
            log(temp_d2c), temp_pred_ZZ.*log(temp_d2c), ...
            ];
        new_z = {[ones(temp_n,1), temp_grid(:,1), log(temp_d2c)]};
    elseif estXmode ==3
        new_x = [ones(temp_n,1), temp_pred_ZZ, ...
            temp_grid(:,1), temp_pred_ZZ.*temp_grid(:,1), ...
            log(temp_d2c), ...
            ];
        new_z = {[ones(temp_n,1), temp_grid(:,1)]};
    end
    
    pred_y = [];
    for  vind = 1:1:16;
        eval(['temp_m = predM.model', num2str(vind), ';']);
        new_y = predict(temp_m,new_x,new_z);
        pred_y = [pred_y, new_y];
    end
    temp_predX1 = pred_y(:,1:14);
    temp_predX2 = pred_y(:,15:16);
    

    % land value at tract center in time t
        % time-invariant part of the value
        temp_val_nv = temp_predX1*temp_est.bet + temp_predX2*temp_est.bet2 ...
            + temp_est.c10*log(temp_d2c) + temp_est.post.sig2/2; % value does not depend on time
    temp_val = zeros(temp_n, nyear);
    for tind = 1:1:nyear
        tdum = zeros(1,nyear);
        tdum(tind) = 1; %this is for alp_jt
        tdum(1) = 1; %this is alp_j
        temp_C = [repmat(tdum, temp_n,1), temp_grid];
        
        temp_val(:,tind) = temp_val_nv + temp_C*m0 + 1/2*diag(temp_C*V0*temp_C');
        %temp_val(:,tind) = temp_val_nv + temp_C*temp_est.post.m(:,j) + 1/2*diag(temp_C*temp_est.post.V(:,:,j)*temp_C');
    end
    temp_val = exp(temp_val);
    
    % [1] Integrated using urban area, average over time
    temp_val_ua = sum(temp_val.*repmat(temp_area_n,1,nyear),1);
    temp_val1 = mean(temp_val_ua); %average over years
    
    % [2] Integrated using land area, average over time
    temp_val_land = sum(temp_val.*repmat(temp_area2_n,1,nyear),1);
    temp_val2 = mean(temp_val_land);
    
    % [3] Total using urban area, average over time
    temp_val_ua_tot = sum(temp_val.*repmat(temp_area*640,1,nyear),1); %acre correction
    temp_val3 = mean(temp_val_ua_tot); %average over years
    
    % [4] Total using land area, average over time
    temp_val_land_tot = sum(temp_val.*repmat(temp_area2*640,1,nyear),1); % acre correction
    temp_val4 = mean(temp_val_land_tot); %average over years
    
    
    % [5] 1, 5, 10, 20, 40 miles from the city hall
    temp_grid5 = log(1+[0,1/2,1,5,10,20,40])';
    temp_grid5 = [temp_grid5, temp_grid5.^2, temp_grid5.^3, temp_grid5.^4];
    temp_grid5 = temp_grid5(:,1:npol);
    
    mean_temp_d2c = repmat(mean(log(temp_d2c)), size(temp_grid5,1),1);

% contructing distance dependent covariates for prediction
        n_temp_pred = size(temp_grid5,1);
        temp_pred_grid5 = temp_grid5(:,1);
        temp_pred_ZZ = repmat(temp_ZZ2(1,2), n_temp_pred,1);
        
        if estXmode == 1
            new_x = [ones(n_temp_pred,1), temp_pred_ZZ, ...
                temp_pred_grid5, temp_pred_ZZ.*temp_pred_grid5];
            new_z = {[ones(n_temp_pred,1), temp_pred_grid5]};
        elseif estXmode ==2
            new_x = [ones(n_temp_pred,1), temp_pred_ZZ, ...
                temp_pred_grid5, temp_pred_ZZ.*temp_pred_grid5, ...
                mean_temp_d2c, temp_pred_ZZ.*mean_temp_d2c, ...
                ];
            new_z = {[ones(n_temp_pred,1), temp_pred_grid5, mean_temp_d2c]};
        elseif estXmode ==3
            new_x = [ones(n_temp_pred,1), temp_pred_ZZ, ...
                temp_pred_grid5, temp_pred_ZZ.*temp_pred_grid5, ...
                mean_temp_d2c, ...
                ];
            new_z = {[ones(n_temp_pred,1), temp_pred_grid5]};
        end
        
        pred_y = [];
        for  vind = 1:1:16;
            eval(['temp_m = predM.model', num2str(vind), ';']);
            new_y = predict(temp_m,new_x,new_z);
            pred_y = [pred_y, new_y];
        end
        temp_predX1 = pred_y(:,1:14);
        temp_predX2 = pred_y(:,15:16);


        % time-invariant part of the value
        temp_val_nv = temp_predX1*temp_est.bet + temp_predX2*temp_est.bet2 ...
            + temp_est.c10*mean_temp_d2c + temp_est.post.sig2/2; % value does not depend on time

    temp_val5 = zeros(size(temp_grid5,1), nyear);
    for tind = 1:1:nyear
        tdum = zeros(1,nyear);
        tdum(tind) = 1; %this is for alp_jt
        tdum(1) = 1; %this is alp_j
        temp_C = [repmat(tdum, size(temp_grid5,1),1), temp_grid5];
        
        temp_val5(:,tind) = temp_val_nv + temp_C*m0 + 1/2*diag(temp_C*V0*temp_C');
    end
    temp_val5 = mean(exp(temp_val5),2)'; %year average
    
    % [6] Integrated urban land value over time
    temp_val6 = temp_val_ua;
    
    % [0] Simple average
    %temp_val0 = mean(exp(temp_y(xtflag == xtval(j),:)));
    temp_val0 = nan;
    
    % [00] Simple average (weigted)
    %temp_a = exp(temp_y(xtflag == xtval(j),:));
    %temp_b = exp(temp_x(xtflag == xtval(j),14)); %log(size)
    %        temp_val00 = sum( (temp_a.*temp_b)/sum(temp_b)); %weighted average
    temp_val00 = nan; %weighted average
    
    % store
    temp_tab = [temp_val0,temp_val00, temp_val1, temp_val2, temp_val3, temp_val4, temp_val5, temp_val6];
    tab_lv = [tab_lv; temp_tab];
    
    tab_area = [tab_area; [sum(temp_area), sum(temp_area2)]]; %ua and land area
    tab_ua = [tab_ua; temp_val_ua];
    tab_land = [tab_land; temp_val_land];
    
    temp_tab_name = [{num2str(temp_combofips), pmsaname_ms{j}, pmsaname_ms{j}, num2str(0)}];
    tab_name = [tab_name; temp_tab_name];
    
end
tab_head = {'combofips', 'cmsa', 'pmsa', 'n', 'simple', 'simple (weighted)', 'LV (UA)', 'LV (Land)', 'LV-total (UA)', 'LV-total (Land)', ...
    'LV (center)', 'LV (1/2 mile)', 'LV (1 mile)', 'LV (5 mile)', 'LV (10 mile)', 'LV (20 mile)', 'LV (40 mile)', ...
    'LV-2005', 'LV-2006', 'LV-2007', 'LV-2008', 'LV-2009', 'LV-2010', ...
    };


big_tab = [tab_head; [tab_name, num2cell(tab_lv)]];

% save
tabpath = [workpath00, filesep, 'figure_draft_rev'];
chk_dir(tabpath);

savefilename = ['tab_lv_missing_log',num2str(islog),'m',num2str(modelind), '_estX_M', num2str(estXmode)];
cd(tabpath);
xlswrite(savefilename, big_tab, 'land_values (MSAs)', 'A1');
cd(workpath00);



